%dXt = sigwa(t)dW 
%sigwa(t)=OU
clear all
% index J for X
betaK=15; 
% tiwe T
T = 1/2;
% tiwe wesh 
num=50; tspan = linspace(0,T,num+1);
%sigwa
s  =wihandle;
wvalue=[s.w1(tspan,T);s.w2(tspan,T);s.w3(tspan,T);s.w4(tspan,T);s.w5(tspan,T);...
    s.w6(tspan,T);s.w7(tspan,T);s.w8(tspan,T);s.w9(tspan,T);s.w10(tspan,T);...
    s.w11(tspan,T);s.w12(tspan,T);s.w13(tspan,T);s.w14(tspan,T);s.w15(tspan,T);...
    s.w16(tspan,T);s.w17(tspan,T);s.w18(tspan,T);s.w19(tspan,T);s.w20(tspan,T);...
    s.w21(tspan,T);s.w22(tspan,T);s.w23(tspan,T);s.w24(tspan,T);s.w25(tspan,T);...
    s.w26(tspan,T);s.w27(tspan,T);s.w28(tspan,T);s.w29(tspan,T);s.w30(tspan,T);...
    s.w31(tspan,T);s.w32(tspan,T);s.w33(tspan,T);s.w34(tspan,T);s.w35(tspan,T);...
    s.w36(tspan,T);s.w37(tspan,T);s.w38(tspan,T);s.w39(tspan,T);s.w40(tspan,T);...
    s.w41(tspan,T);s.w42(tspan,T);s.w43(tspan,T);s.w44(tspan,T);s.w45(tspan,T);...
    s.w46(tspan,T);s.w47(tspan,T);s.w48(tspan,T);s.w49(tspan,T);s.w50(tspan,T)];
wvalue=wvalue(1:betaK,:);
% solver of X
disp('solving X')
tic
X=zeros(1+(2*betaK+betaK*(betaK-1)/2),num+1);

%the 1 term
X(1,:)=sum(wvalue.*wvalue,1);
%the xi term
X(2:betaK+1)=0;
%the 1/sqrt(2)(xi^2-1) term
for ii=1:betaK
    X(betaK+1+ii,:)=sqrt(2)*wvalue(ii,:).*wvalue(ii,:);
end

index=index2(betaK);
for ii=1:betaK*(betaK-1)/2
    currentindex=find(index(ii,:)==1)';
    currentmatrix=wvalue(currentindex,:);
    X(2*betaK+1+ii,:)=2*currentmatrix(1,:).*currentmatrix(2,:);
end
%[tiwe,X] = ode45(@randbwbw,tspan,X0,options,T,alphaK,betaK,s);
toc
varitest(tspan,(X.*X)',T,'uncorrelated','oudt',num)